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We introduce an efficient Langevin metfiod to study biiinear Fermionic Hamiitonians interacting 
witfi ciassicai fieids. Our metfiod is suitabfe for very farge systems and offers fiigfi accuracy To 
demonstrate tfie metfiod, we study compfex non-copfanar cfiiraf spin textures on tfie triangufar 
Kondo fattice modef. We afso expfore non-equifibrium mesoscafe pfiysics sucfi as cfiiraf domain 
coarsening and Z2 vortex annifiifation. 



Lattice modefs of Fermions interacting witfi cfassicaf 
fiefds encompass a wide range of pfiysics. Exampfes in 
condensed matter incfude Kondo fattice (KL) modefs [Ij, 
Faficov-Kimbaff modefs [2j and Bogofiubov-De Gennes 
equations [3j. Tfiese modefs are numericaffy cfiaffenging 
because eacfi triaf cfiange to tfie cfassicaf fiefd introduces 
a new singfe particfe Fermion matrix to be diagonafized. 

KL modefs in particufar fiave recentfy emerged as a 
ricfi source of pfiysics. Tfiey are effective for d- and /- 
efectron compounds and describe tfie coupfing between 
itinerant efectrons and effectivefy cfassicaf focafized mag- 
netic moments. Tfie cofossaf magnetoresistance (CMR) 
effect [H [5] spurred tfie study of KL modefs at farge 
coupfing, wfiere a ferromagnetic transition of tfie focaf- 
ized magnetic moments feads to a dramatic cfiange in 
resistivity. Severaf fast Monte- Car fo (MC) metfiods fiave 
been devefoped and appfied to tfiis ferromagnetic transi- 
tion [6 - 9j , wfiere onfy moderate temperature and numer- 
icaf precision is needed. 

Recent interest in tfie KL modef fias sfiifted to excit- 
ing non-copfanar spin textures of tfie focaf magnetic fiefd. 
For exampfe, Sl^yrmion fattices fiave recentfy been ob- 
served HT] witfi spatiaffy modufated spin textures on 
fengtfi scafes up to 0.1/im. Anotfier exampfe is cfiiraf 
spin textures, in wfiicfi efectrons interact witfi a fiuge ef- 
fective magnetic fiefd (~ lO^T) due to tfie Berry pfiase. 
Tfiis anomafous Haff effect fias been predicted in tfie KL 
modef witfi triangufar fattice [12j and fias been experi- 
mentaffy observed in, e.g., tfie compounds Pr2lr207 [T3] 
and UCu5 [14J. Unfil^e CMR, tfie numericaf study of 
tfiese exotic pfiases typicaffy requires very farge system 
sizes and fiigfi numericaf precision. 

In tfiis paper we introduce a Langevin sampfing 
metfiod tfiat is very efficient — tfie cost scafes linearly 
witfi system size — at fiigfi accuracy. Our metfiod is based 
on a gradient transformation [15] of tfie l^ernef pofynomiaf 
metfiod (KPM) [16]. Befow we outfine our metfiod and 
tfien use it to study equifibrium cfiiraf spin textures in 
tfie triangufar KL modef. Our fattices are farge enougfi 
to uncover interesting non-equifibrium effects sucfi as cfii- 
raf domain coarsening and Z2 vortex dynamics. In tfiis 
way, we bridge tfie gap between quantum atomic scafe 
and mesoscafe pfiysics. 

We present our metfiod in tfie generaf context of a 
bifinear Fermionic Hamiftonian coupfed to continuous. 



cfassicaf degrees of freedom ^, 

n = clAij{(j))cj, 



(1) 



witfi sparse matrix A. We worl^ at fixed temperature 
and cfiemicaf potentiaf ja. Tfie partition function 
is a trace over cfassicaf and Fermionic degrees of free- 
dom, Z = TT(f)TTcex.-p[—/3{l-L — fi^.clci)]. Evafuating 
tfie Fermionic trace yiefds Z = Tr^ exp(— wfiere 



e)de 
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is tfie effective (free) energy of configuration (/), p{e) = 
J(e — ejy{(j))) is tfie density of states of and 
/(e) = —f3~^ log{l + exp[— /3(e — fi)]}. Tfie energy con- 
tains effective fong-range many-body interactions of tfie 
cfassicaf fiefd. 

A l^ey difficufty in MC sampfing tfie cfassicaf fiefd is tfie 
cafcufation of AF in response to cfianges in 0. KPM es- 
timates tfie density of states p{e) using a series of Cfieby- 
sfiev pofynomiafs Tm{^) truncated at order M ^T6\ \T7\ . 
Tfie cost of KPM is finear in system size N for sparse ma- 
trix A. We state directfy tfie recursive KPM procedure 
to construct an unbiased stocfiastic estimate of F 
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{r m = 

Ar m = l (5) 

2Aam-i - OLm-2 m> 1 

Here, r is a random cofumn vector wfiose components 
satisfy (r*rj) = 6ij. We draw from a uniform distribu- 
tion |r^p = 1 to guarantee tfiat tfie Cfiebysfiev moments 
Pm are exact wfien A is diagonaf. Tfie coefficients 

Cm = J - e2) (2 - Jo,m) gruTm (e) /(e)de 

are independent of A. Tfie Jacl^son l^ernef 
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is chosen to damp Gibbs oscillations yet retain high ac- 
curacy. KPM requires that the eigenvalues of A have 
magnitude less than 1, which can usually be achieved 
without loss of generality as a rescaling of energy. 

We emphasize that there are two sources of error in 
KPM: (1) truncation at finite order M and (2) stochas- 
tic estimation by averaging over finitely many random 
vectors r. Both are well controlled and will be discussed 
below. 

To sample fields {(/)} from the Boltzmann distribution, 
P[(j)] (X exp{—/3F [(/)]), we apply the overdamped Langevin 
equation. In discretized form, 

dF , 

(i>i(t + At) - (i>i(t) = -At-- + y^2(3-^Atr]i{t), (6) 

where r]i{t) are uncorrelated Gaussian random variables 
with unit variance and t is a fictitious time. The Langevin 
approach simultaneously updates all components Ef- 
ficient and accurate estimation of the gradient V iF = 
dF/d(j)i is crucial. We exclude inertial terms from the 
Langevin equation because they would amplify errors in 
the estimate of VF. 

The technique of automatic differentiation with "re- 
verse accumulation" [15j ensures that, by careful appli- 
cation of the chain rule, we can transform the KPM pro- 
cedure to estimate F, Eqs. [3}|5) into one that estimates 
VF at the same cost. We perform this transformation 
analytically and state only the final result. 



— -— = /3o;iao;j + 2 Pm-i(^m;j (7) 

J m=l 

The row vectors are given by reverse recursion, from 
m = M — 2 down to m = 0, 

/3m-2 = (8) 

Cm+ir^ + 2l3m+lA - /3^+2 (9) 

The desired gradient is dF/d(j)i = 
^j^i{dF/dAi^i){dA}^i/d(pi). The sequence of vectors 
am are the same as in the original KPM, but are here 
required in reverse order. We recalculate them as needed 
using am = 2Aam-\-i — The recursion begins with, 

aM-i and aM-2, which are available at the end of the 
original KPM procedure. 

Note that the procedure to estimate V^F, Eqs. [7]-[8j 
has computational cost equivalent to the original KPM 
procedure to calculate F. 

The gradient calculation also inherits the approxima- 
tion errors of KPM, controlled by two parameters: the 
truncation order M of the Chebyshev series, and the dy- 
namical stochastic error which is the integration time- 
step At per KPM random vector. The KPM estimated 
density of states p{e) is resolved to order Ae/M, where 
Ae = Cmax — Cmin is the spau of extremal eigenvalues. The 



parameter M should be chosen large enough to resolve 
the physically relevant features of the density of states. 
Finite stochastic error z > acts much like an additional 
noise term in the Langevin dynamics, effectively rescaling 
its magnitude by an amount T ^ T + ATgff . By scaling, 
we expect ATeff ~ 2:, independent of system size. How- 
ever, if A = Aq is diagonal then there is no stochastic er- 
ror and ATeff = 0. For perturbations A = Aq^JAi with 
J <C 1 the stochastic error remains small; by symmetry, 
the lowest order term is ATeff ^ J'^z. The parameter z 
should be chosen small enough that ATeff <^ T for the 
smallest relevant temperature scale T. 

Our Langevin sampling remains efficient at high accu- 
racy: the cost to integrate the Langevin equation one unit 
of time is 0{NM / z). To compare to Metropolis MC with 
local updates, we assume that one unit of Langevin inte- 
gration time roughly corresponds to a full MC sweep in 
which all N lattice sites are visited. Trial changes A(j)i are 
accepted with a probability that depends on the change 
in energy, AT. Brute force exact diagonalization of ma- 
trix A requires 0{N^) operations, and a full MC sweep 
costs 0{N^). This may be reduced to 0{N^) by track- 
ing the response of the spectrum to low-rank changes in 
A [8l[T8]. Further acceleration is possible with KPM ap- 
proximation. A non-stochastic Green's function method 
reduces the cost of a fuh MC sweep to 0{MN'^) [9j. In 
an alternative approach, if A contains only local coupling 
in d dimensions, a full MC sweep may be performed at 
cost C)(M^+i7V) [7J. In the simulat ions below we will use 
d = 2, N = 100^ M 1000, and z 0.02, for which 
our Langevin approach outperforms existing methods by 
2 orders of magnitude or more. Another linear cost al- 
gorithm is based on hybrid Monte-Carlo [6j, but a direct 
comparison is difficult because the cost of precision at 
low temperatures is unclear. 

We apply our method with the triangular KL model, 
whose Hamiltonian is 

H = - ^ tijc\^Cj^ -J^Sj- c]^a^^Cj^, (10) 

where c^^(cjo-) is the creation (annihilation) operator of 
an electron with spin a on site j, Sj is a classical Heisen- 
berg spin with |S^ | = 1, and o-^^ = (cr^^, cr^^, cr^^) is 
a vector of Pauli matrices. The hopping coefficients are 
tij = tex.-p{—iOij) when i and j are nearest neighbor lat- 
tice sites, and t^j = otherwise. In the following, we fix 
the energy scale by taking t ^ 1, and the spatial scale by 
taking the lattice spacing to 1. A possibly non-uniform 
phase Oij = BzZ • x Xj/2 models orbital coupling to 
the dimensionless magnetic field B = BzZ^ where is 
the position of lattice site i. 

Martin and Batista argue, by perfect nesting of the 
Fermi surface, that the chiral 3g configuration (Fig. ^) 
is the ground state at 3/4 electron filling fraction with 
small coupling [12j. This state is of special interest, as it 
exhibits a spontaneous quantum Hall effect. Variational 
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Figure 1. Three competing periodic 2x2 spin-textures in the 
triangular Kondo lattice model at 3/4 electron filling fraction. 
The (a) 1^, (b) 2^, and (c) 3^ ("all-out") phases are named 
according to their number of reciprocal lattice vectors. The 
3q phase maximizes chirality x = Sa x Sb • Sc = ±4/3'^^^ of 
triangular plaquettes and gives rise to an spontaneous quan- 
tum Hall effect at 1/4 and 3/4 fillings. The Iq and 2q phases 
break rotational symmetry of the triangular lattice. 



calculation on the 2x2 plaquette also predicts stabil- 
ity of the 3q state [19j. However, unconstrained Monte- 
Carlo study of this phase has not yet been achieved due 
to severe numerical difficulties. Very low temperatures, 
T < 0.001 and small couplings J < 0.3 are required 
to stabilize 3^^. The resolution of the numerical method 
must be very high to resolve the small gap in the density 
of states (of width ~ J^). Furthermore, the 3q state is 
stabilized by a susceptibility that diverges like log^ so 
very large lattices sizes are required {N ^ 100^). These 
challenges imply that this system provides a rigorous test 
of our Langevin method. 

We choose J 0.2 and /i = 1.947, and TV = 
100^ There are three competing phases, 3q^ 2q^ 
and Iq (Fig. [T]) with very similar energy densities, 
-4.15552,-4.15550,-4.15525, respectively. We use 
KPM based Langevin sampling with M = 1000 and z = 
0.02. At M = 1000, KPM estimates of energy differences 
are accurate to order 10~^. With z = 0.02 the effective 
Langevin temperature is increased by ATeff « 0.0002. 

In Fig. [2|a) we observe melting of the chiral 3q ground 
state. The chirality x abruptly disappears at a first or- 
der phase transition at T ^ 0.0010. This transition, 
however, is not to a paramagnetic phase. To distinguish 
the phases, 3g, 2g, and Iq we consider the (unordered) 
set of nearest-neighbor spin-spin correlations C = {(Sx • 

Sx+(i,o)), {^^■K+{i,V3)/2)^ {^^-K+{-i,V3)/2)} averaged 
over lattice sites x. Pure 3q^ 2g, and Iq phases would 
yield Csq = {-1/3,-1/3,-1/3}, C2, = {0,0,1}, and 
Ciq = {-1,-1,1}. The latter two states have broken 
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Figure 2. Phase diagram of the triangular Kondo lattice 
model at Hund coupling J = 0.2 and chemical potential 
/i = 1.947, corresponding to filling fraction ^ 3/4. (a) At 
low temperatures the preferred 3^ phase is identified by its 
non-zero chirality x — (Sa x St • Sc). Very large system sizes 
are required to stabilize the 3^ phase; N — 100^ (circles) is 
sufficient, but N — 60^ (squares) is not. (b) Spin-spin correla- 
tion functions (Sa-Sb) for three nearest-neighbor orientations. 
Three first order phase transitions are apparent: (1) 3^ to 2^ 
at T 0.0010, (2) 2q to 1^ at T = 0.0017, and (3) Iq to 
paramagnet at T = 0.0029. 



bond symmetry (specifically, the 3-fold rotational sym- 
metry of the triangular lattice). 

Fig. |2|b) plots the three elements of C as a function 
of temperature. At T = we find C = Csq as expected. 
We now observe three first order transitions at tempera- 
tures T = 0.0010, 0.0017, and T = 0.0029 to the 2q, Iq, 
and paramagnetic phases, respectively. The 2q phase is 
identified by its correlation set C, which has two zero el- 
ements and one negative element. In the Iq phase, C has 
one positive element and two negative (symmetric) ele- 
ments. To avoid equilibration issues in the above data, 
we used initial conditions with explicitly broken chiral 
symmetry, x > 0? 

We now investigate the dynamical, non-equilibrium 
process by which 3q chiral symmetry breaking occurs at 
low temperatures. We use our Langevin dynamics to 
study the phase ordering kinetic of chiral domains fol- 
lowing a quench from infinite to zero temperature. 

Langevin "time", properly speaking, is fictitious. How- 
ever, in the spirit of time dependent Ginzburg-Landau 
(TDGL), we expect energetic relaxation to capture the 
qualitative, large-scale aspects of phase ordering [20j. A 
point of comparison is Model A dynamics in the Ho- 
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Figure 3. Phase ordering in the 200 x 200 triangular Kondo lattice model following a quench from infinite to zero temperature. 
The color gradient, ranging from red to blue, is the local chirality. Langevin times are measured in units of r = 1.28 x 10^. (a) 
J = 3, /i = —3.2 1/4 filling). Domain coarsening with strong anisotropy is observed. Z2 vortices appearing as white dots 
rapidly annihilate each other, (b) J = 0.2, /a = 1.947 (^3/4 filling). The system is dynamically trapped in a complex, robust 
metastable state, (c) J = 0.2^ /i = 1.947, Bz = Stt/VSN 3/4 filling). An external field breaks chiral symmetry and the 
system rapidly evolves to the 3^ ground state. A Z2 vortex is identified by the winding of a Burger's circuit (green) in SO (3) 
space, a filled projective sphere in the axis-angle representation. 



henberg and Halperin classification [21j, the prototypi- 
cal TDGL model for phase-ordering of a non-conserved 
scalar order parameter. One expects Model A to apply 
generally to this "kinetic Ising" universality class. In- 
terestingly, we find that the Kondo lattice model has ef- 
fective long-range many-body interactions that introduce 
new dynamical features not present Model A. 

First we consider the case of ~ 1/4 filling with J = 3 
and /i = — 3.2. Previous work found a robust 3q phase at 
system sizes up to = 16^ [22j. We use our Langevin 
dynamics to study the ordering dynamics ai N = 200^ 
with accuracy parameters M = 500, z = 0.005. Fig- 
ure [2|a) shows the evolution of chirality ranging from 
red (positive) to blue (negative). The coarsening of chiral 
domains is analogous to Model A, but we observe strong 
anisotropy of domain walls. The dynamics slows at large 
times, consistent with a characteristic length scale that 
grows as £ - t^/^ [23j. 

Next we consider ~ 3/4 filling, with J = 0.2 and 
jii = 1.947. We use accuracy parameters M = 500 



and z = 0.02. The ordering dynamics is displayed in 
Fig. [2|b). A new dynamical feature appears: the chi- 
ral domains evolve into a remarkable pattern which is a 
very long lived metastable state. This is a reproducible 
phenomenon. At higher temperatures this metastable 
pattern could be annealed to the pure 3q phase, but in 
experimental practice it is easier to explicitly break chiral 
symmetry with an applied external field [13]. With 
the smallest possible orbital coupling for a periodic lat- 
tice, Bz = 87r/200v/3, we find that system rapidly reaches 
the uniform chiral 3^^ phase, shown in Fig. [2|c). 

Topological defects, visible as small white dots, are ap- 
parent at both 1/4 and 3/4 filling. These defects are Z2 
vortices associated with winding of the SO {3) topologi- 
cal manifold. A Z2 vortex is enlarged in Fig.[2|c), second 
and third panels. One of the four 3q spin sub-lattices are 
shown. We can understand this defect by constructing 
a closed Burger's circuit (draw in green) that encircles 
it. Every point on the green circuit is identified with an 
element in 50(3), plotted as a trajectory in the fourth 
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panel. The SO {3) manifold, in the axis-angle represen- 
tation, is a filled-sphere with antipodal points identified. 
This Burger's circuit has winding number 1 because it 
wraps SO {3). The vortex is Z2 because the only homo- 
topically distinct winding numbers are and 1. Conse- 
quently, any pair of Z2 vortices may annihilate. In the 
ordering dynamics of Fig.[2|a), we observe many vortices 
annihilating with each other and with domain walls. 

In conclusion, we have introduced a numerical method 
to study the broad class of Hamiltonians that couple 
Fermions to classical degrees of freedom. Our method 
is highly accurate and efficient, enabling the study of 
complex systems at unprecedented size. Large system 
sizes may be necessary to resolve logarithmic divergences 
associated with nesting of the Fermi surface. In the tri- 
angular Kondo lattice model at 3/4 filling, we found that 
lattices of size N = 100^ with six digits of precision are 
required. Large system sizes also allow us to bridge the 
gap between quantum and mesoscopic physics. With sys- 
tems of size N = 200^ we are able to probe chiral domain 
dynamics, metastable trapping, and Z2 vortex dynamics, 
effects typically inaccessible in standard approaches. 
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